%--------------------------------------------------------------------------
% computes the "neutral stability" curves that are predicted by the method
% of multiple scales
%
% input arguments:
% kk is the range of k values to plot
% tt are the times at which to construct the curves
% fig is a flag for plotting the results
%--------------------------------------------------------------------------


function [kk, Ma, M, small_k] = exact_neutral_stab(kk, tt, figs)
if figs == 1
    close all
end

p = params;

if (nargin == 1)
    tt = [0, 1, 10, 100];
end

Lambda_c = zeros(length(tt), length(kk));
Ma_c = Lambda_c;

for i = 1:length(tt)
    t = tt(i);
    h_bar = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * t) / (p.beta - 1))))) * (p.beta - 1);
    
    h = h_bar;
    v_0 = (1 - p.beta) * (1 - 1 / h_bar);    

    c_bar = p.beta + v_0 - p.delta / 2 * (1 - p.beta) * (p.beta + v_0);
    gam = p.delta * (1 - 2 * c_bar);
    
    delta = p.delta;
    beta = p.beta;
    
    for j = 1:length(kk)
        k = kk(j);

        Ma(i,j) = (0.12e2 * (-1 + beta) * (beta + v_0) * k * (0.5e1 / 0.2e1 + k ^ 2 * h ^ 2) * h * delta ^ 2 * cosh(k * h) ^ 3 - 0.48e2 * (k ^ 2 * (c_bar - 0.1e1 / 0.2e1) * h ^ 2 + delta * (-1 + beta) * (beta + v_0) / 0.4e1) * sinh(k * h) * delta * cosh(k * h) ^ 2 - 0.12e2 * k * ((k ^ 2 * delta ^ 2 * (-1 + beta) * (beta + v_0) * h ^ 2 - 0.2e1 * k ^ 2 * h + delta ^ 2 * (-1 + beta) * (beta + v_0)) * sinh(k * h) ^ 2 + 0.5e1 / 0.6e1 * (k ^ 2 * (delta * (-1 + beta) * (beta + v_0) + 0.12e2 / 0.5e1 - 0.24e2 / 0.5e1 * c_bar) * h ^ 2 + 0.9e1 / 0.5e1 * delta * (-1 + beta) * (beta + v_0)) * delta) * h * cosh(k * h) - 0.4e1 * (k ^ 2 * delta ^ 2 * (-1 + beta) * (beta + v_0) * h ^ 2 + 0.6e1 * k ^ 2 * h + 0.3e1 * delta ^ 2 * (-1 + beta) * (beta + v_0)) * sinh(k * h) * k ^ 2 * h ^ 2) / delta / (-1 + beta) / (beta + v_0) / (0.6e1 * cosh(k * h) * sinh(k * h) ^ 2 * k ^ 3 * h ^ 3 - 0.3e1 * sinh(k * h) + 0.3e1 * cosh(k * h) ^ 2 * sinh(k * h) + 0.3e1 * h * cosh(k * h) * k + 0.3e1 * k ^ 2 * h ^ 2 * sinh(k * h) + 0.3e1 * k ^ 3 * h ^ 3 * cosh(k * h) - 0.6e1 * cosh(k * h) ^ 3 * k ^ 3 * h ^ 3 - 0.3e1 * cosh(k * h) ^ 3 * h * k + 0.2e1 * k ^ 4 * h ^ 4 * sinh(k * h));
        small_k(i,j) = -(2 * (-11 * delta * beta - 11 * delta * v_0 + 11 * delta * beta ^ 2 + 11 * delta * beta * v_0 - 80 * c_bar + 40) / (beta + v_0) / (-1 + beta) / h ^ 2 / k ^ 2) - 0.4e1 / 0.63e2 / h * (-85 * h * delta ^ 2 * beta - 85 * h * delta ^ 2 * v_0 + 85 * h * delta ^ 2 * beta ^ 2 + 85 * h * delta ^ 2 * beta * v_0 - 904 * h * delta * c_bar + 452 * h * delta + 1260) / (beta + v_0) / delta / (-1 + beta);
      
        M(i,j) = -(-(-k * sinh(k * h) - delta * (0.1e1 - (2 * beta) - 0.2e1 * delta * beta * (1 - beta) * ((h - 0.1e1) / beta / h / delta - (h - 0.1e1 + beta) / beta / h / 0.2e1 + 0.1e1 / 0.6e1)) * cosh(k * h)) / (-h / 0.12e2 + 0.1e1 / 0.12e2 - beta / 0.12e2) * beta * h ^ 3 * k ^ 2 / (0.6e1 * h * cosh(k * h) ^ 2 * k - k ^ 3 * h ^ 3 - 0.3e1 * cosh(k * h) * sinh(k * h) - 0.3e1 * k * h - 0.3e1 * k ^ 2 * h ^ 2 * cosh(k * h) * sinh(k * h)) * (-k * h + cosh(k * h) * sinh(k * h)) - 0.2e1 * delta ^ 2 * (beta ^ 2) * cosh(k * h) + 0.2e1 * delta ^ 2 * beta * cosh(k * h) + 0.2e1 * k * delta ^ 2 * beta * sinh(k * h) * h - 0.2e1 * k * delta ^ 2 * (beta ^ 2) * sinh(k * h) * h) / sinh(k * h) / k / h;

M(i,j) =(0.2e1 * (-k * sinh(k * h) - delta * (0.1e1 - (2 * beta) - 0.2e1 * delta * beta * (1 - beta) * ((h - 0.1e1) / beta / h / delta - (h - 0.1e1 + beta) / beta / h / 0.2e1 + 0.1e1 / 0.6e1)) * cosh(k * h)) / (-h / 0.12e2 + 0.1e1 / 0.12e2 - beta / 0.12e2) * beta * h ^ 3 * k ^ 2 * (-k * h + cosh(k * h) * sinh(k * h)) + 0.18e2 * h * cosh(k * h) * k * delta ^ 2 * beta - 0.18e2 * h * cosh(k * h) * k * delta ^ 2 * (beta ^ 2) + 0.12e2 * cosh(k * h) * sinh(k * h) ^ 2 * k * delta ^ 2 * beta * h - 0.12e2 * cosh(k * h) * sinh(k * h) ^ 2 * k * delta ^ 2 * (beta ^ 2) * h + 0.12e2 * k ^ 3 * h ^ 3 * cosh(k * h) * sinh(k * h) ^ 2 * delta ^ 2 * beta - 0.12e2 * k ^ 3 * h ^ 3 * cosh(k * h) * sinh(k * h) ^ 2 * delta ^ 2 * (beta ^ 2) - 0.10e2 * k ^ 3 * h ^ 3 * delta ^ 2 * (beta ^ 2) * cosh(k * h) + 0.10e2 * k ^ 3 * h ^ 3 * delta ^ 2 * beta * cosh(k * h) + 0.4e1 * k ^ 4 * h ^ 4 * delta ^ 2 * beta * sinh(k * h) - 0.4e1 * k ^ 4 * h ^ 4 * delta ^ 2 * (beta ^ 2) * sinh(k * h) + 0.12e2 * k ^ 2 * h ^ 2 * delta ^ 2 * beta * sinh(k * h) - 0.12e2 * k ^ 2 * h ^ 2 * delta ^ 2 * (beta ^ 2) * sinh(k * h) - 0.12e2 * cosh(k * h) ^ 2 * sinh(k * h) * delta ^ 2 * (beta ^ 2) + 0.12e2 * cosh(k * h) ^ 2 * sinh(k * h) * delta ^ 2 * beta - 0.30e2 * cosh(k * h) ^ 3 * h * k * delta ^ 2 * beta + 0.30e2 * cosh(k * h) ^ 3 * h * k * delta ^ 2 * (beta ^ 2) - 0.12e2 * k ^ 3 * cosh(k * h) ^ 3 * h ^ 3 * delta ^ 2 * beta + 0.12e2 * k ^ 3 * cosh(k * h) ^ 3 * h ^ 3 * delta ^ 2 * (beta ^ 2)) / (-0.6e1 * k ^ 3 * h ^ 3 * cosh(k * h) * sinh(k * h) ^ 2 - 0.3e1 * cosh(k * h) ^ 2 * sinh(k * h) + 0.3e1 * sinh(k * h) - 0.3e1 * h * cosh(k * h) * k - 0.2e1 * k ^ 4 * h ^ 4 * sinh(k * h) - 0.3e1 * k ^ 2 * h ^ 2 * sinh(k * h) - 0.3e1 * k ^ 3 * h ^ 3 * cosh(k * h) + 0.3e1 * cosh(k * h) ^ 3 * h * k + 0.6e1 * k ^ 3 * cosh(k * h) ^ 3 * h ^ 3);
% M(i,j) = (0.2e1 * (-k * sinh(k * h) - delta * (0.1e1 - (2 * beta) - 0.2e1 * delta * (1 - beta) * (h - 0.1e1) / h) * cosh(k * h)) / (-h / 0.12e2 + 0.1e1 / 0.12e2 - beta / 0.12e2) * beta * h ^ 3 * k ^ 2 * (-k * h + cosh(k * h) * sinh(k * h)) + 0.18e2 * h * cosh(k * h) * k * delta ^ 2 * beta - 0.18e2 * h * cosh(k * h) * k * delta ^ 2 * (beta ^ 2) + 0.12e2 * cosh(k * h) * sinh(k * h) ^ 2 * k * delta ^ 2 * beta * h - 0.12e2 * cosh(k * h) * sinh(k * h) ^ 2 * k * delta ^ 2 * (beta ^ 2) * h + 0.12e2 * k ^ 3 * h ^ 3 * cosh(k * h) * sinh(k * h) ^ 2 * delta ^ 2 * beta - 0.12e2 * k ^ 3 * h ^ 3 * cosh(k * h) * sinh(k * h) ^ 2 * delta ^ 2 * (beta ^ 2) - 0.10e2 * k ^ 3 * h ^ 3 * delta ^ 2 * (beta ^ 2) * cosh(k * h) + 0.10e2 * k ^ 3 * h ^ 3 * delta ^ 2 * beta * cosh(k * h) + 0.4e1 * k ^ 4 * h ^ 4 * delta ^ 2 * beta * sinh(k * h) - 0.4e1 * k ^ 4 * h ^ 4 * delta ^ 2 * (beta ^ 2) * sinh(k * h) + 0.12e2 * k ^ 2 * h ^ 2 * delta ^ 2 * beta * sinh(k * h) - 0.12e2 * k ^ 2 * h ^ 2 * delta ^ 2 * (beta ^ 2) * sinh(k * h) - 0.12e2 * cosh(k * h) ^ 2 * sinh(k * h) * delta ^ 2 * (beta ^ 2) + 0.12e2 * cosh(k * h) ^ 2 * sinh(k * h) * delta ^ 2 * beta - 0.30e2 * cosh(k * h) ^ 3 * h * k * delta ^ 2 * beta + 0.30e2 * cosh(k * h) ^ 3 * h * k * delta ^ 2 * (beta ^ 2) - 0.12e2 * k ^ 3 * cosh(k * h) ^ 3 * h ^ 3 * delta ^ 2 * beta + 0.12e2 * k ^ 3 * cosh(k * h) ^ 3 * h ^ 3 * delta ^ 2 * (beta ^ 2)) / (-0.3e1 * cosh(k * h) ^ 2 * sinh(k * h) - 0.6e1 * k ^ 3 * h ^ 3 * cosh(k * h) * sinh(k * h) ^ 2 + 0.3e1 * sinh(k * h) - 0.3e1 * h * cosh(k * h) * k - 0.2e1 * k ^ 4 * h ^ 4 * sinh(k * h) - 0.3e1 * k ^ 2 * h ^ 2 * sinh(k * h) - 0.3e1 * k ^ 3 * h ^ 3 * cosh(k * h) + 0.3e1 * cosh(k * h) ^ 3 * h * k + 0.6e1 * k ^ 3 * cosh(k * h) ^ 3 * h ^ 3);

    end
end


if figs
    plot(kk, p.delta * p.beta * (1 - p.beta) * Ma);
    ylim([0, max(get(gca, 'ylim'))]);
end



% h_bar = fzero(@(h_bar) fun(h_bar, k, p), 0.9);
% tau = 1 - h_bar + (p.beta - 1) * log10((p.beta - 1 + h_bar) / p.beta);
% t = tau / p.delta
% 
% function v = fun(h_bar, k, p)
% v_0 = (1 - p.beta) * (1 - 1 / h_bar);
% A = w_coeffs(1, h_bar, k, p) / p.Ma;
% Lambda = 0.12e2 * (p.delta * (1 - 2 * p.beta) * cosh(k * h_bar) + sinh(k * h_bar) * k) * k ^ 3 / ((0.2e1 * A(1) * h_bar ^ 3 * k ^ 3 - 0.6e1 * k * h_bar * A(1) + 0.3e1 * A(3) * h_bar ^ 2 * k ^ 2) * cosh(k * h_bar) ^ 2 + 0.6e1 * sinh(k * h_bar) * (A(1) * h_bar ^ 2 * k ^ 2 - A(3) * h_bar * k + A(1)) * cosh(k * h_bar) - 0.2e1 * sinh(k * h_bar) ^ 2 * (A(1) * h_bar ^ 3 * k ^ 3 + 0.3e1 * k * h_bar * A(1) - 0.3e1 / 0.2e1 * A(3) * h_bar ^ 2 * k ^ 2 - 0.3e1 / 0.2e1 * A(3)));
% 
% v = p.Ma - Lambda * h_bar / p.delta / (1 - p.beta - v_0) / (p.beta + v_0);